% 猪肝样本的物理属性
rho = 1050; % 猪肝密度 (kg/m^3)
c = 3500; % 比热容 (J/kg·K)
kappa = 0.5; % 热导率 (W/m·K)

% 电极的加热区域体积 (椭球体积)
a = 0.01; % 长轴 (m)
b = 0.008; % 短轴 (m)
c_depth = 0.03; % 深度 (m)
V = (4/3) * pi * a * b * c_depth; % 椭球体积 (m^3)
P = 30; % 功率 (W)

% 离散化网格
r_max = 0.05; % 最大r
z_max = 0.05; % 最大z
dr = 0.0001; % r的步长
dz = 0.0001; % z的步长
dt = 0.01; % 时间步长
t_max = 15; % 最大时间

% 网格尺寸
r = -r_max:dr:r_max;
z = -z_max:dz:z_max;
t = 0:dt:t_max;

% 初始化温度场
T = zeros(length(r), length(z), length(t));

% 初始条件
T(:,:,1) = 25; % 初始温度 25℃

% 修正热源（椭球区域）的影响
% 电极加热区域椭球形状：r, z 只在椭球区域内加热
for n = 1:length(t)-1
    for i = 2:length(r)-1
        for j = 2:length(z)-1
            % 电极加热区域，假设电极加热区域在椭圆内
            if (r(i)^2/a^2 + z(j)^2/b^2) <= 1  % 椭圆区域内
                % 离散化热传导方程
                T(i,j,n+1) = T(i,j,n) + (kappa*dt/(rho*c)) * ((T(i+1,j,n) - 2*T(i,j,n) + T(i-1,j,n))/dr^2 + (T(i,j+1,n) - 2*T(i,j,n) + T(i,j-1,n))/dz^2) + (P/V) * dt;
            else
                % 在椭圆外的区域，温度不变（或按边界条件设置）
                T(i,j,n+1) = T(i,j,n);  % 也可以设置为一个常温或边界温度
            end
        end
    end
end

% 选择某个时间步进行展示
function plotTemperatureField(T, r, z, time_step)
    % 获取对应时间步的温度场
    T_at_time = T(:,:,time_step);
    T_at_time1 = T_at_time / max(T_at_time(:)) * 100;  % 归一化温度场

    % 校正r和z坐标
    r_corrected = r * (2.27/0.0075) * (2.7/2);  
    z_corrected = z * (2.27/0.0075) * (2.7/2);

    % 绘制温度场
    figure;
    imagesc(r_corrected, z_corrected, T_at_time1);  % 绘制二维温度场
    axis xy;  % 设置坐标轴方向正确
    xlabel('r (mm)');
    ylabel('z (mm)');
    title(['Temperature Distribution at t = ', num2str(time_step), ' s']);
    colorbar;  % 显示颜色条

    % 移除z轴，隐藏z轴标签
    zlim([min(T_at_time(:)), max(T_at_time(:))]);  % 设置z轴的范围
    set(gca, 'ztick', []);  % 隐藏z轴刻度
    set(gca, 'zticklabel', []);  % 隐藏z轴刻度标签

    axis equal;  % 保证坐标轴比例相同
end
function plotIntensity(T, r, z, time_step)
    % 获取对应时间步的温度场
    T_at_time = T(:,:,time_step);
    T_at_time1 = T_at_time / max(T_at_time(:)) * 100;  % 归一化温度场

    % 获取第500行的数据（r方向上的强度）
    intensity_data = T_at_time1(500, :);

    % 绘制强度图
    figure;
    plot(r, intensity_data, 'LineWidth', 2);
    xlabel('r (mm)');
    ylabel('Temperature (%)');
    title(['Temperature Intensity at Row 500 (t = ', num2str(time_step), ' s)']);
    grid on;
end

% time_step = 300;
% T_at_time=T(:,:,time_step)
% T_at_time1 = T_at_time / max(T_at_time(:)) * 100;
% % 绘制最终的温度场
% % 使用 imagesc 绘制颜色映射图
% figure;
% r_corrected = r * (2.27/0.0075)* (2.7/2);  % 横坐标 (r) 进行校正
% z_corrected = z * (2.27/0.0075)* (2.7/2);;  % 纵坐标 (z) 进行校正
% % 
% % imagesc(r_corrected, z_corrected, new_image);  % 绘制二维温度场
% imagesc(r_corrected, z_corrected, T_at_time1);  % 绘制二维温度场
% axis xy;  % 设置坐标轴方向正确
% xlabel('r (mm)');
% ylabel('z (mm)');
% title(['Temperature Distribution at t = ', num2str(time_step), ' s']);
% colorbar;  % 显示颜色条
% 
% 
% % 移除 z 轴，隐藏 z 轴标签
% zlim([min(T_at_time(:)), max(T_at_time(:))]);  % 设置 z 轴的范围
% set(gca, 'ztick', []);  % 隐藏 z 轴刻度
% set(gca, 'zticklabel', []);  % 隐藏 z 轴刻度标签
% 
% axis equal;  % 保证坐标轴比例相同
% new_image = zeros(1001 + 228, 1001);
% % 归一化每一层温度
% % 归一化
% time_step = 900;
% T_at_time=T(:,:,time_step)
% T_at_time1 = T_at_time / max(T_at_time(:)) * 100;
% 
% % 将 T_at_time 叠加228次，每次向上平移一行
% for i = 1:228
%     new_image(i:i+size(T_at_time1, 1)-1, :) = new_image(i:i+size(T_at_time1, 1)-1, :) + T_at_time1;
% 
%     % 如果超过150，就设置为150
% 
% end
% new_image = new_image / max(new_image(:)) * 100;
% % 显示最终结果
% 
% % 绘制最终的温度场
% figure;
% r_corrected = r * (2.27/0.0075)* (2.7/2);  % 横坐标 (r) 进行校正
% z_corrected = z * (2.27/0.0075)* (2.7/2);;  % 纵坐标 (z) 进行校正
% 
% imagesc(r_corrected, z_corrected, new_image);  % 绘制二维温度场
% caxis([30 40]);
% 
% axis xy;  % 设置坐标轴方向正确
% xlabel('r (mm)');
% ylabel('z (mm)');
% title(['Temperature Distribution at t = ', num2str(time_step*0.01), ' s']);
% colorbar;  % 显示颜色条
% 
% % 设置 z 轴的范围
% zlim([0, 100]);  % 归一化后最大温度为100
% set(gca, 'ztick', []);  % 隐藏 z 轴刻度
% set(gca, 'zticklabel', []);  % 隐藏 z 轴刻度标签
% 
% axis equal;  % 保证坐标轴比例相同



% % 绘制最终的温度场
% % 使用 imagesc 绘制颜色映射图
% figure;
% imagesc(r, z, T_at_time);  % 绘制二维温度场
% axis xy;  % 设置坐标轴方向正确
% xlabel('r (m)');
% ylabel('z (m)');
% title(['Temperature Distribution at t = ', num2str(time_step), ' s']);
% colorbar;  % 显示颜色条
% 
% 
% % 移除 z 轴，隐藏 z 轴标签
% zlim([min(T_at_time(:)), max(T_at_time(:))]);  % 设置 z 轴的范围
% set(gca, 'ztick', []);  % 隐藏 z 轴刻度
% set(gca, 'zticklabel', []);  % 隐藏 z 轴刻度标签
% 
% axis equal;  % 保证坐标轴比例相同
